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ABSTRACT 

We use a set of high-resolution cosmological N-body simulations to investigate the inner 
mass profile of galaxy-sized cold dark matter (CDM) halos. These simulations extend the thor- 
ough numerical convergence study presented in Paper I of this series (Power et al. 2003), and 
demonstrate that the mass profile of CDM halos can be robustly estimated beyond a minimum 
converged radius of order r conv ~ 1 /i _1 kpc in our highest resolution runs. The density profiles 
of simulated halos become progressively shallow from the virial radius inwards, and show no sign 
of approaching a well-defined power-law behaviour near the centre. At r con v, the logarithmic 
slope of the density profile is steeper than the asymptotic p oc r" 1 expected from the formula 
proposed by Navarro, Frenk, and White (1996), but significantly shallower than the steeply di- 
vergent p oc r -1 - 5 cusp proposed by Moore et al. (1999). We perform a direct comparison of the 
spherically-averaged dark matter circular velocity (V c ) profiles with rotation curves of low surface 
brightness (LSB) galaxies from the samples of de Blok et al. (2001b), de Blok and Bosma (2002), 
and Swaters et al. (2003). Most (about two-thirds) LSB galaxies in this dataset are roughly 
consistent with CDM halo V c profiles. However, about one third of LSBs in these samples feature 
a sharp transition between the rising and flat part of the rotation curve that is not seen in the 
V c profiles of CDM halos. This discrepancy has been interpreted as excluding the presence of 
cusps, but we argue that it might simply reflect the difference between circular velocity and gas 
rotation speed likely to arise in gaseous disks embedded within realistic, triaxial CDM halos. 

matter - galaxies: formation - galaxies: spiral - 

on the scale free nature of the gravitational accre- 
tion process and suggested that halo density pro- 
files might be simple power laws (Gunn and Gott 
1972; Fillmore and Goldreich 1984; Hoffman and 
Shaham 1985; White and Zaritsky 1992). Cosmo- 
logical N-body simulations, however, failed to con- 
firm these analytic expectations. Although power- 
laws with slopes close to those motivated by the 
theory were able to describe some parts of the halo 
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galaxies: kinematics and dynamics 

1. Introduction 

The structure of dark matter halos and its re- 
lation to the cosmological context of their forma- 
tion has been studied extensively over the past 
few decades. Early analytic calculations focused 
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density profiles, even early simulations found that 
significant deviations from a single power-law be- 
haviour were present in most cases (Frenk et al. 
1985, 1988; Quinn et al. 1986; Dubinski and Carl- 
berg 1991; Crone et al. 1994). Further simulation 
work, indeed, concluded that power-law fits were 
inappropriate, and that, properly scaled, dark ha- 
los spanning a wide range in mass and size are well 
fit by a "universal" density profile (Navarro et al. 
1995, 1996, 1997, hereafter NFW): 

PNFwM = , -, u f S -j—r^- (1) 

(r/r s )(l + r/r s ) 2 

One characteristic of this fitting formula is that 
the logarithmic slope, (3(r) = —d\og p/d\ogr = 
(1 + 3y)/(l + y) (where y — r/r s is the radius in 
units of a characteristic scale radius, r s ), increases 
monotonically from the centre outwards. The den- 
sity profile steepens with increasing radius; it is 
shallower than isothermal inside r s , and steeper 
than isothermal for r > r s . Another important 
feature illustrated by this fitting formula is that 
the profiles are "cuspy" (/3o = (3{r = 0) > 0): the 
dark matter density (but not the potential) di- 
verges formally at the centre. 

Subsequent work has generally confirmed these 
trends, but has also highlighted potentially impor- 
tant deviations from the NFW fitting formula. In 
particular, Fukushige and Makino (1997, 2001), 
as well as Moore and collaborators (Moore et al. 
1998, 1999; Ghigna et al. 2000), have reported 
that NFW fits to their simulated halos (which 
had much higher mass and spatial resolution than 
the original NFW work) underestimate the dark 
matter density in the innermost regions (r < r s ). 
These authors proposed that the disagreement was 
indicative of inner density "cusps" steeper than 
the NFW profile and advocated a simple modifica- 
tion to the NFW formula so that /?o = 1-5 (rather 
than 1.0). 

The actual value of the asymptotic slope, /3o, 
is still being hotly debated in the literature (Jing 
et al. 1995; Klypin et al. 2001; Taylor and Navarro 
2001; Navarro 2003; Power et al. 2003; Fukushige 
et al. 2003), but there is general consensus that 
CDM halos are indeed "cuspy". This has been 
recognized as an important result, since the rota- 
tion curves of many disk galaxies, and in particular 
of low surface brightness (LSB) systems, appear to 
indicate the presence of an extended region of con- 



stant dark matter density: a dark matter "core" 
(Flores and Primack 1994; Moore 1994; Burkert 
1995; Blais-Ouellette et al. 2001; de Blok et al. 
2001a,b). 

Unfortunately, rotation curve constraints are 
strongest just where numerical simulations are 
least reliable. Resolving CDM halos down to 
the kpc scales probed by the innermost points of 
observed rotation curves requires extremely high 
mass and force resolutions, as well as careful inte- 
gration of particle orbits in the central, high den- 
sity regions of halos. This poses a significant com- 
putational challenge that has been met in very few 
of the simulations published to date. 

This difficulty has meant that rotation curves 
have usually been compared with extrapolations 
of the simulation data that rely heavily on the ap- 
plicability and accuracy of fitting formulae such as 
the NFW profile to regions that may be severely 
compromised by numerical artifact. This practice 
does not allow for halo-to- halo variations, or for 
temporary departures from the "average" profile 
to be properly taken into account when trying to 
model the observational datasets. 

Finally, the theoretical debate on the asymp- 
totic central slope of the dark matter density pro- 
file, /3q, has led at times to unwarranted empha- 
sis on the very inner region of the rotation curve 
datasets, rather than on an proper appraisal of 
the data over its full radial extent. For example, 
de Blok et al. (2001a, b) attempt to derive con- 
straints on 00 from the innermost few points of 
their rotation curves, and conclude that /?o ~ 
for most galaxies in their sample. However, their 
analysis focuses on the regions most severely af- 
fected by non-circular motions, seeing, misalign- 
ments and slit offsets, and other effects that limit 
the accuracy of circular velocity estimates based 
on long-slit spectra. It is perhaps not surprising, 
then, that other studies have disputed the conclu- 
siveness of these findings. 

Indeed, an independent analysis of data of sim- 
ilar quality by Swaters et al. (2003) (see also van 
den Bosch et al. 2000) concludes that "cuspy" dark 
matter halos with (3o <^ 1 are actually consistent 
with their data. This disagreement is compounded 
by the results of the latest cosmological N-body 
simulations (Power et al. 2003, hereafter P03), 
which find scant evidence for a well defined value 
of Po m simulated CDM halos. Given these dif- 
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Acuities, focusing the theoretical or observational 
analysis on /3o seems unwise. 

In this paper, we improve upon previous work 
by comparing directly the actual results of the sim- 
ulations with the full radial extent of the rotation 
curves of LSB galaxies. We present results from a 
set of seven galaxy-sized dark matter halos, each 
of which has been carefully simulated at various 
resolution levels in order to ascertain the numer- 
ical convergence of our results. This allows us 
to test rigorously the P03 convergence criteria, as 
well as to assess the cusp-core discrepancy through 
direct comparison between observation and simu- 
lations. A companion paper (Navarro et al. 2003) 
addresses the issue of universality of CDM halo 
structure using simulations that span a wide range 
of scales, from dwarf galaxies to galaxy clusters. 

The outline of this paper is as follows. In §2 
we introduce our set of simulations and summa- 
rize briefly our numerical methods. The seven 
galaxy-sized halos that form the core of our sam- 
ple have been simulated at various resolutions, and 
we use them in §3 to investigate the robustness of 
the P03 numerical convergence criteria. The den- 
sity profiles of these halos are presented and com- 
pared with previous work in §4. In §5 we compare 
the halo V c profiles with the LSB rotation curve 
datasets of de Blok et al. (2001b), de Blok and 
Bosma (2002), and Swaters et al. (2003). Our 
main conclusions and plans for future work are 
summarized in §6. 

2. The Numerical Simulations 

We have focused our analysis on seven galaxy- 
sized dark matter halos selected at random 
from two different cosmological N-body simula- 
tions of large periodic boxes with comoving size 
i box = 32.5/i _1 Mpc and 35.325 /i _1 Mpc, respec- 
tively. Each of these "parent" simulations has 
Nbox — 128 3 particles, and adopts the currently 
favoured flat, low-density "concordance" ACDM 
world model, with Qq = 0.3, Qa = 0.7, and ei- 
ther h = 0.65 (runs labelled Gl, G2 and G3) or 
h = 0.7 (G4, G5, G6, and G7, see Table 1). The 
power spectrum in both simulations is normalized 
so that the linear rms amplitude of fluctuations 
on spheres of radius 8 /i _1 Mpc is as — 0.9. 

All halos (Gl to G7) have been re-simulated 
at three or four different mass resolution levels; 



each level increases the number of particles in the 
halo by a factor of 8, so that the mass per particle 
has been varied by a factor 512 in runs G1-G3, 
and by a factor 64 in runs G4-G7 (see Table 1). 
All of these runs focus numerical resources on the 
Lagrangian region from where each system draws 
its mass, whilst approximating the tidal field of 
the whole box by combining distant particles into 
groups of particles whose mass increases with dis- 
tance from the halo. This resimulation technique 
follows closely that described in detail in P03 and 
in Navarro et al. (2003), where the reader is re- 
ferred to for full details. For completeness, wc 
present here a brief account of the procedure. 

Halos selected for resimulation are identified at 
z = from the full list of halos with circular veloc- 
ities in the range (150, 250) km s _1 in the parent 
simulations. All particles within a sphere of radius 
3 ?*200 2 centred in each halo are then traced back 
to the initial redshift configuration (z.i = 49). 

The region defined by these particles is typi- 
cally fully contained within a box of size L s ^ ox ~ 
5 /i _1 Mpc, which is loaded with A sbox = 32 3 , 64 3 , 
128 3 , or 256 3 particles. Particles in this new high- 
resolution region are perturbed with the same 
waves as in the parent simulation as well as with 
additional smaller scale waves up to the Nyquist 
frequency of the high resolution particle grid. Par- 
ticles which do not end up within 3 r2oo of the se- 
lected halo at z = are replaced by lower resolu- 
tion particles which replicate the tidal field acting 
on the high resolution particles. This resampling 
includes some particles within the boundaries of 
the high resolution box, and therefore the high res- 
olution region defines an asymmetrical "amoeba- 
shaped" three-dimensional volume surrounded by 
tidal particles whose mass increases with distance 
from this region. 

A summary of the numerical parameters and 
halo properties is given in Table 1. This table also 
includes reference to 12 further runs, four of them 
corresponding to dwarf galaxy halos and eight of 
them to galaxy cluster-sized halos. These systems 
have been simulated only at the highest resolution 
(-Wsbox = 256 3 ), and therefore are not included 

2 We define the "virial radius", r200, as the radius of a sphere 
of mean density 200 times the critical value for closure, 
Pcrit = 3 H 2 /8ttG, where H is Hubble's constant. We pa- 
rameterize the present value of Hubble's constant H by 
H = 100 h km s" 1 kpc" 1 
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in our convergence analysis. These runs are dis- 
cussed in detail in a companion paper (Navarro 
et al. 2003). 

Some simulations were performed with a fixed 
number of timesteps for all particles using Stadel 
and Quinn's parallel N-body code PKDGRAV (Stadel 
2001), while others used the N-body code GADGET 
(Springel et al. 2001). The GADGET runs allowed 
for individual timesteps for each particle assigned 
using either the RhoSgAcc or EpsAcc criterion (see 
P03 for full details). The halo labelled Gl in this 
paper is the same one selected for the numerical 
convergence study presented in P03. Although 
PKDGRAV also has individual timestepping capabili- 
ties, we have chosen not to take advantage of these 
for the simulations presented in this paper. We 
note that P03 finds only a modest computational 
gain due to multi-stepping schemes provided that 
the softening parameter is properly chosen. 

The softening parameter (fixed in comoving co- 
ordinates) for each simulation (with the exception 
of Gl/256 3 , see P03) was chosen to match the "op- 
timal" softening suggested by P03: 



e opt 



4 r 20 o 



N. 



1/2 



(2) 



where N200 is the number of particles within r2oo 
at z = 0. This softening choice minimizes the 
number of timesteps required for convergence re- 
sults by minimizing discreteness effects in the force 
calculations whilst ensuring adequate force resolu- 
tion. 

At z = 0, the mass within the virial radius, 
M200, of our galaxy-sized halos ranges from ~ 
10 12 h- x M Q to - 3 x 10 12 h~ l M Q , correspond- 
ing to circular velocities, V 2 oo = (GM 2 ao/r 2 oa) 1/2 , 
in the range 160 kms" 1 to 230 kms -1 . 

3. Numerical Convergence 
3.1. Criteria 

P03 propose three different conditions that 
should be satisfied in order to ensure conver- 
gence in the circular velocity profile. According to 
these criteria, convergence to better than 10% in 
the spherically-averaged circular velocity, V c (r), 
is achieved at radii which satisfy the following 
conditions: 



greater than the size of the timestep At: 

tcirc(f200) V ^0 J 

where t denotes the age of the universe, 
which is by definition of the order of the 
circular orbit timescale at the virial radius, 

£circ(f20C>)- 

2. Accelerations do not exceed a characteristic 
acceleration, a e , determined by V200 and the 
softening length e: 



a(r) = 



V c 2 (r) 



<,a t = 0.5 



V 2 

v 200 



(4) 



where a(r) is the mean radial accelera- 
tion experienced by particles at a distance 
r from the centre of the system, a{r) = 
GM{r)/r 2 = V c 2 (r)/r. 

3. Enough particles are enclosed such that the 
local collisional relaxation timescale t Te \ ax (r) 
is longer than the age of the universe 3 : 



^rclax 



(r) 



/200 N(r) (p(r) 



icirc(r-20o) 8 lniV(r) \p cr it 



-1/2 



£1 

(5) 

where N(r) is the number of particles and 
p(r) is the mean density within radius r. 

For "optimal" choices of the softening and 
timestep, as well as for the typical number of 
particles in our runs, we find that criterion (3) 
above is the strictest one. 

The number of high-resolution particles thus ef- 
fectively defines the "predicted" converged radius, 
r conv , beyond which, according to P03, circular ve- 
locities should be accurate to better than 10%. We 
emphasize that this accuracy criterion applies to 
the cumulative mass profile; convergence in prop- 
erties such as local density estimates, p(r), typ- 
ically extends to radii significantly smaller than 



3 We adopt a slightly more conservative criterion than P03, 
who require t ro i ax > 0.6 t c irc(»*20()). 



1. The local orbital timescale t c [ rc (r) is much 
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3.2. Validating the Convergence Criteria 

We assess the validity of the convergence crite- 
ria listed above by comparing the mass profile of 
the highest resolution run corresponding to each 
halo with those obtained at lower resolution. Fig- 
ure 1 illustrates the procedure. From top to bot- 
tom, the three panels in this figure show, as a 
function of radius, the circular orbit timescale, 
the mean radial acceleration, and the relaxation 
timescale, respectively, for the four runs corre- 
sponding to halo G3. The small arrows at the 
bottom of each panel indicate the choice of gravi- 
tational softening for each run. The dotted curves 
in the top and middle panel show the best fit NFW 
profile to the converged region of the highest res- 
olution 7V s box = 256 3 run. 

The "converged radius" corresponding to each 
criterion is determined by the intersection of the 
horizontal dashed lines in each panel with the 
"true" profile, which we shall take to be that of 
the highest resolution run (shown in solid black in 
Figure 1). Clearly, the strictest criterion is that 
imposed by the relaxation timescale (the dotted 
vertical lines in the lower panel show the converged 
radius corresponding to this criterion). This sug- 
gests, for example, that the lowest-resolution G3 
run (with A^box = 32 3 , shown in solid blue), 
should start to deviate from the converged profiles 
roughly at r <~ 0.1 r 2 oo- Indeed, this appears to 
be the radius at which this profile starts to "peel 
off" from the highest resolution one, as shown in 
the top two panels of Figure 1. Increasing the 
number of high resolution particles by a factor of 
eight typically brings the converged radius inwards 
by a factor of <~ 2.4. For the medium-resolution 
run (A^ s box = 64 3 , shown in solid green), r conv is 
predicted to be ~ 0.04 r 2 oo, which, again, coin- 
cides well with the radius inside which departures 
from the converged profile are apparent. Simi- 
larly, r conv <~ 0.017 r 2 oo for the high-resolution 
(-^Vsbox — 128 3 ) run (shown in red). 

The density and circular velocity profiles cor- 
responding to the four G3 runs are shown in Fig- 
ure 2. Panels on the left show the profiles down 
to the radius that contains 50 particles, whereas 
those on the right show the profiles restricted to 
r ^> r conv . Figure 2 illustrates two important re- 
sults alluded to above: (i) both p(r) and V c (r) 
converge well at r ;> r conv , and (ii) convergence in 



p(r) extends to radii smaller than r conv . Indeed, 
the top-left panel shows that our choice of r conv 
is rather conservative when applied to the density 
profile. Typically, densities are estimated to bet- 
ter than 10% down tor~ 0.6r conv . 

How general are these results? Figure 3 com- 
pares the minimum "converged" radius predicted 
by the P03 criteria , r conv , with r w % vc , the actual 
radius where circular velocities in the lower resolu- 
tion runs deviate from convergence by more than 
10%. In essentially all cases, r conv <; r w % vc , indi- 
cating that the P03 criteria are appropriate, albeit 
at times somewhat conservative. We list our r conv 
estimates for all runs in Table 1. 

4. Halo Structure and Fitting Formulae 

The dotted curves in Figure 2 show the best 
NFW fits to the density and circular velocity pro- 
file of the highest resolution run. The dashed lines 
correspond to the best fit adopting the modifica- 
tion to the NFW profile advocated by Moore et al. 
(1999), 

/ \ Pm_ (R , 

PMooie[r) - (r/r M y*(l + (r/r M y*y W 

These fits are obtained by straightforward \ 2 min- 
imization in two parameters, r s or tm, and the 
characteristic density p s or pm- The profiles are 
calculated in bins of equal width in logr, and the 
fits are performed over the radial range r conv < 
t < f 200 • Equal weights are assigned to each radial 
bins because the statistical (Poisson) uncertainty 
in the determination of the mass within each bin 
is negligible (each bin contains thousands of par- 
ticles) so uncertainties are completely dominated 
by systematic errors whose radial dependence is 
difficult to assess quantitatively. 

The best fits to p(r) and V c (r) shown in Fig- 
ure 2 are obtained independently from each other. 
Values of the concentration parameter, cnfw = 
^200 As, for the best fit NFW profiles are 6.4 and 
5.3 for fits to the density and circular velocity 
profile, respectively; the Moore et al concentra- 
tions, CMooro = r20o/fM, are 3.0 and 2.9 for the 
best fits to p(r) and V c (r), respectively. Over the 
converged region, r ;> r conv , both the NFW and 
Moore et al profiles appear to reproduce reason- 
ably well the numerical simulation results. Indeed, 
no profile in the G3 runs deviates by more than 
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10% in V c or 30% in p(r) from the best fits ob- 
tained with either eq. 1 or eq. 6. More substantial 
differences are expected only well inside r con v, but 
these regions are not reliably probed by the simu- 
lations. 

This suggests that either the NFW or Moore et 
al profile may be used to describe the structure of 
ACDM halos outside ~ 1% of the virial radius, but 
also implies that one should be extremely wary of 
extrapolations inside this radius. One intriguing 
feature of Figure 2 is that the Moore et al formula 
appears to fit the G3 density profiles better than 
NFW but that V c profiles are somewhat better ap- 
proximated by NFW (see also P03). This suggests 
that neither formula captures fully and accurately 
the radial dependence of the structure of ACDM 
halos. 

4.1. The radial dependence of the loga- 
rithmic slope 

This view is confirmed by the radial depen- 
dence of the logarithmic slope of the density profile 
(3(r) — —d\ogp/d\ogr, which is shown in the top- 
left panel of Figure 4 for all the high-resolution 
runs, and compared with the predictions of the 
NFW (solid line) and Moore et al (dashed line) 
formulae. 

Logarithmic slopes are calculated by numer- 
ical differentiation of the density profile, com- 
puted in radial bins of equal logarithmic width 
(A logr/r 2 oo — 0.2). The slope profiles in Figure 4 
are normalized to r_2, the radius where (3{r) takes 
the "isothermal" value of 2. 4 In this and all sub- 
sequent figures, profiles are shown only down to 
the minimum converged radius r conv . This corre- 
sponds typically to a radius r conv ~ 0.006 r 2 oo, or 
about 1-2 h~ 1 kpc for halos simulated at highest 
resolution (see Table 1). 

The top left panel of Figure 4 shows that halos 
differ from the NFW and Moore et al formula in 
a number of ways: 

• there is no obvious convergence to an asymp- 
totic value of the logarithmic slope at the 
centre; the profile gets shallower all the way 
down to the innermost radius reliably re- 
solved in our runs, r conv . 

4 r_2 is in this sense equivalent to the scale radius r s of the 
NFW profile. 



• the slope at r conv is significantly shallower 
than the asymptotic value of (3o — 1.5 advo- 
cated by Moore et al. (1999). The shallowest 
value measured at r conv is (3 ~ 1, and the av- 
erage over all seven halos is (3 ~ 1.2. 

• Most halo profiles become shallower with ra- 
dius more gradually than predicted by the 
NFW formula; at r ~ 0.1 r^i the average 
slope is ~ —1.4, whereas NFW would predict 

1.18. The NFW density profile turns 

over too sharply from p cx r~ 3 to p cx 
compared to the simulations. 

In other words, the Moore et al profile appears 
to fit better the inner regions of the density pro- 
file of some ACDM halos (see bottom-left panel of 
Figure 4) not because the inner density cusp di- 
verges as steeply as (3q = 1.5, but rather because 
its logarithmic slope becomes shallower inwards 
less rapidly than NFW. 

It is important to note as well that there is sig- 
nificant scatter from halo to halo, and that two 
of the seven density profiles are actually fit bet- 
ter by the NFW formula. Are these global devi- 
ations from a "universal" profile due to substruc- 
ture? We have addressed this question by remov- 
ing substructure from all halos and then recom- 
puting the slopes. Substructure is removed by 
first computing the local density at the position 
of each particle, pi, using a spline kernel similar 
to that used in Smoothed Particle Hydrodynamics 
(SPH) calculations 5 . Then, we remove all particles 
whose densities are more than 2 standard devia- 
tions above the spherically-averaged mean density 
at its location. (The mean and standard devi- 
ation are computed in bins of equal logarithmic 
width, Alogr/r 2 oo — 0.01). The procedure is it- 
erated until no further particles are removed. The 
remaining particles form a smoothly distributed 
system that appears devoid of substructure on all 
scales. We find that density profiles are smoother 
after the removal of substructure but that most 
of the variation in the overall shapes of the pro- 
files remains. We conclude that the presence of 
substructure is not directly responsible for the ob- 
served scatter in the shape of halo density profiles. 

See http: //www-hpec . astro .Washington . edu/tools/smooth .html 
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4.2. Comparison with Other Work 

Are these conclusions consistent with previous 
work? To explore this issue, we have computed 
the logarithmic slope profile of three CDM halos 
run by Moore and collaborators. The halos we 
have re-analyzed are the Milky Way- and M31- 
like galaxy halos of the Local Group system from 
Moore et al. (1999) and the LORES version of the 
"Virgo" cluster halo from Ghigna et al. (2000). 
The z ~ 0.1 output of the Local Group simulation 
was provided to us by the authors, whereas the 
Virgo cluster was re-run using initial conditions 
available from Moore's website 6 . The Virgo clus- 
ter run used the same N-body code as the origi- 
nal simulation (PKDGRAV) but was run with a fixed 
number of timesteps (12800). A run with 6400 
timesteps was also carried out and no differences 
in the mass profiles were detected. The number 
of particles within the virial radius is 1.2 x 10 6 , 
1.7 x 10 6 , and 5.0 x 10 5 , for the Milky Way (MW), 
M31 and LORES Virgo cluster halos, respectively. 

Figure 4 shows the logarithmic slope (upper 
right panel) and density (lower right panel) pro- 
files corresponding to these halos, plotted down 
to the minimum converged radius r conv . No ma- 
jor differences between these simulations and ours 
are obvious from these panels. It is clear, for 
example, that at the innermost converged point, 
the slope of the density profile of the two Local 
Group halos is significantly shallower than r^ 15 , 
and shows no signs of having converged to a well 
defined power-law behaviour. There is some evi- 
dence for "convergence" to a steep cusp (r~ 1A ) in 
the LORES Virgo cluster simulation but the dy- 
namic range over which this behaviour is observed 
is rather limited. The Virgo cluster run thus ap- 
pears slightly unusual when compared with other 
systems in our ensemble. Although our reanaly- 
sis confirms the conclusion of Moore et al. (1998, 
1999) that this system appears to have a steeply 
divergent core, this does not seem to be a general 
feature of ACDM halos. 

Our results thus lend support to the conclu- 
sions of Klypin et al. (2001), who argues that there 
is substantial scatter in the inner profiles of cold 
dark matter halos. Some are best described by the 



http://www.nbody.net. We note that all of these runs were 
evolved in an f2o = 1 cosmogony, rather than the ACDM 
scenario wc adopt in this paper. 



NFW profile whereas others are better fit by the 
Moore et al formula, implying that studies based 
on a single halo might reach significantly biased 
conclusions. 

Finally, we note that deviations from either fit- 
ting formula in the radial range resolved by the 
simulations, although significant, are small. Best 
NFW/Moore et al fits are typically accurate to 
better than ~ 20% in circular velocity and ~ 40% 
in density, respectively. We discuss in a compan- 
ion paper (Navarro et al. 2003) the constraints 
placed by our simulations on extrapolations of 
these formulae to the inner regions as well as on 
the true asymptotic inner slope of ACDM halo 
density profiles. 

5. Halo Circular Velocity Profiles and LSB 
Rotation Curves 

As discussed in § 1, an important discrepancy 
between the structure of CDM halos and the mass 
distribution in disk galaxies inferred from rotation 
curves has been noted repeatedly in the literature 
over the past decade (Moore 1994; Flores and Pri- 
mack 1994; Burkert 1995; McGaugh and de Blok 
1998; Moore et al. 1999; van den Bosch et al. 2000; 
Cote ct al. 2000; Blais-Ouellette et al. 2001; van 
den Bosch and Swaters 2001; Jimenez et al. 2003). 
In particular, the shape of the rotation curves of 
low surface brightness (LSB) galaxies has been 
identified as especially difficult to reconcile with 
the "cuspy" density profiles of CDM halos. 

Given the small contribution of the baryonic 
component to the mass budget in these galax- 
ies, the rotation curves of LSB disks are expected 
to trace rather cleanly the dark matter potential, 
making them ideal probes of the inner structure of 
dark matter halos in LSBs. Many of these galax- 
ies are better fit by circular velocity curves arising 
from density profiles with a well defined constant 
density "core" rather than the cuspy ones inferred 
from simulations, a result that has prompted calls 
for a radical revision of the CDM paradigm on 
small scales (see e.g., Spergel and Steinhardt 2000) 

It is important, however, to note a number of 
caveats that apply to the LSB rotation curve prob- 
lem. 

• Many of the early rotation curves where 
the disagreement was noted were unduly af- 
fected by beam smearing in the HI data 
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(Swaters et al. 2000). For example, van den 
Bosch et al. (2000) argue that, once beam 
smearing is taken into account, essentially all 
HI LSB rotation curves are consistent with 
cuspy halo profiles. The observational situa- 
tion has now improved substantially thanks 
to higher-resolution rotation curves obtained 
from long-slit H a observations (see, e.g., Mc- 
Gaugh et al. 2001; de Blok et al. 2001a; Swa- 
ters et al. 2000, 2003). We shall restrict our 
analysis to these rotation curves in what fol- 
lows. 

• Strictly speaking, the observational dis- 
agreement is with the fitting formulae, 
rather than with the actual structure of sim- 
ulated CDM halos. As noted in the previous 
section, there are systematic differences be- 
tween them, so it is important to confirm 
that the disagreement persists when LSB 
rotation curves are contrasted directly with 
simulations. 

• Finally, it must be emphasized that the ro- 
tation curve problem arises when compar- 
ing rotation speeds of gaseous disks to the 
spherically-averaged circular velocity pro- 
files of dark matter halos. Given that CDM 
halos are expected to be significantly non- 
spherical (Barnes and Efstathiou 1987; War- 
ren et al. 1992; Jing et al. 1995; Thomas 
ct al. 1998; Jing and Suto 2002), some dif- 
ferences between the two are to be expected. 
It is therefore important to use the full 3D 
structure of CDM halos to make predictions 
regarding the rotation curves of gaseous 
disks that may be compared directly to ob- 
servation. We shall neglect this complex 
issue in this paper, but plan to explore in 
detail the rotation curves of gaseous disks 
embedded in such asymmetric potentials in 
future papers of this series. 

We may avoid many of these uncertainties by 
comparing directly the circular velocity profiles of 
our simulated halos with the observational data. 
This procedure has the advantage of retaining the 
diversity in the shapes of halo profiles that is of- 
ten lost when adopting a simple analytic fitting 
formula. In addition, we consider circular velocity 
profiles only down to the innermost converged ra- 



dius, thereby eliminating uncertainties about the 
reliability of the profile at very small radii. 

We begin the analysis by emphasizing the im- 
portance of taking into account the changes in 
the central halo mass profile induced by accre- 
tion events. Indeed, these may trigger and sustain 
departures from the "average" profile that may 
be detectable in the rotation curves of embedded 
gaseous disks. We shall then describe a simple 
characterization of rotation curve shapes that may 
be applied to both observational and simulation 
data. This enables a direct and quantitative as- 
sessment of the "cusp" versus "core" problem as 
it applies to the most recent LSB datasets. 

5.1. Evolution of the Inner Mass Profile 

Systematic — and at times substantial — changes 
in the inner circular velocity profile are induced by 
accretion events during the assembly of the halo, 
even when such these events might contribute only 
a small fraction of mass to the inner regions. These 
transients may increase substantially the scatter in 
the shape of the V c profiles and they ought to be 
taken into account when comparing with observa- 
tion. 

This is illustrated in Figure 5, which shows the 
evolution of the mass and circular velocity profiles 
of halo Gl. The top panel of this figure shows 
the evolution of the mass enclosed within 8, 10, 
20 kpc (physical), and r 2 oo, as a function of red- 
shift. Although the mass inside 20 kpc increases 
by less than ~ 25% since z = 1, there are signif- 
icant (<~ 50 — 60%) fluctuations during this time 
caused by the tidal effects of orbiting substruc- 
ture and accretion events. Most noticeable is a 
major merger at z ~ 0.7, which affects the mass 
profile down to the innermost reliably resolved ra- 
dius, ~ 2 kpc. 

The effect of these fluctuations on the circular 
velocity profile is shown in the bottom panel of 
Figure 5. Here we show the inner 20 kpc of the 
circular velocity profile before (z = 1.1), during 
(z = 0.48) and after [z = 0) a major accretion 
event. Substantial changes in the shape of the V c 
profile are evident as the halo responds to the in- 
falling substructure. Note that the changes persist 
over timescales of order ;> 1 Gyr, exceeding the 
circular orbital period at r = 2, 10, and 20 kpc 
(~ 0.13, 0.34, and 0.58 Gyr, respectively). These 
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relatively long lasting changes thus would likely be 
reflected in the dynamics of a disk present at the 
centre of the halo. 

5.2. LSB Rotation Curves 

Could the evolutionary effects discussed above 
be responsible, at least in part, for the constant 
density cores inferred from the rotation curves of 
LSB and dwarf galaxies? Since it is nearly impos- 
sible to tailor a simulation to reproduce individ- 
ual galaxies in detail, it is important to adopt a 
simple characterization of the rotation curves that 
allows for a statistical assessment of the disagree- 
ment between halo V c profiles and observation. 
We have thus adopted a three-parameter fitting 
formula commonly used in observational work to 
describe optical rotation curves (Courteau 1997): 



Here Vb is a velocity scale, x = r t /r, where r t 
is a scale radius, and the dimensionless parame- 
ter 7 describes the overall shape of the curve. The 
larger the value of 7 the sharper the turnover from 
the "rising" to the "flat" region of the velocity 
curve. Eq. 7 is flexible enough to accommodate 
the shape of essentially all rotation curves in the 
samples we consider here. We note that this for- 
mula has three 7 free parameters, one more than 
the NFW profile. 

We have applied this fitting formula to the 
HI/Ha rotation curve datascts of de Blok et al. 
(2001b) 8 , de Blok and Bosma (2002) 9 , and Swa- 
ters et al. (2003) 10 , hereafter B01, B02, and S03, 
respectively. Fits to the rotation curves and V c 
profiles are obtained through straightforward x 2 - 
minimization, adopting the error estimates pro- 
vided by the authors. 

The B01 sample consists of 26 LSB galaxies, 
the B02 sample consists of 24 LSB galaxies, and 
the S03 sample contains 10 dwarf galaxies and 5 

7 An additional factor of (1 + x) 13 was used by Courteau 
(1997) to improve fits to ~ 10% of the galaxies in his sample 
that exhibit a drop-off in the outer part of the rotation 
curve. For simplicity, we have not included this parameter 
in our fits. 

8 http : //www . astro . umd . edu/^ssm/data 

9 f tp : //cdsarc . u-strasbg . f r/cats/ J/A+A/385/816 

°http : //www . robswork . net /data 



LSB galaxies. The smoothed rotation curves of 
B01 were derived by folding approaching and re- 
ceding velocities about the centre of the galaxy 
(defined by the peak in the continuum emission) 
and using a spline-fitting procedure followed by 
rebinning to a bin width of 2". Error estimates 
were calculated as the quadratic sum of an ob- 
servational error component caused by measure- 
ment uncertainties in the raw data points, and an 
additional error component due to differences be- 
tween the approaching and receding velocities in 
the bin, as well as noncircular motions (defined 
as the difference between the mean velocity and 
the velocity of the spline fit). The rotation curve 
data provided by B02 are unsmoothed, with er- 
rors which include measurement, inclination, and 
asymmetry uncertainties. The rotation curves of 
S03 were derived by averaging receding and ap- 
proaching velocity data in 2" bins; error estimates 
were defined by half the quadratic sum of the aver- 
age error of points in the bin and half the difference 
between the maximum and minimum velocities in 
the bin. 

The fit parameters for the combined set of 65 
rotation curves are given in Table 2. In this ta- 
ble, r max refers to the radius where the rotation 
curve reaches the maximum, V max . The outermost 
radius with reliable data is listed as r oute r- Fig- 
ure 6 shows a selected sample of rotation curves, 
together with the best fits obtained with eq. 7. 
The top, middle, and bottom rows include galaxies 
from the S03, B02, and B01 samples, respectively, 
arranged from left to right in order of increasing 
value of the shape parameter 7. 

It is important to note that the rotation curves 
of B02 and S03 differ significantly from those of 
B01. The smoothing applied to the B01 dataset is 
clear in this figure, especially when compared with 
the less-processed B02 and S03 datasets. Because 
the B01 curves have been smoothed, the individ- 
ual values quoted for the rotation speed are corre- 
lated, and the error estimates lack clear statistical 
meaning. This is confirmed by the extremely low 
values of the formal reduced \ 2 obtained for the 
best fits with eq. 7 (see Table 2 and Figure 7): 16 
out of 26 B01 galaxies have x 2 ed < 0.1 (and 5 have 
Xred < 0-01), whereas all galaxies in S03's sample 
have Xred > 0.1. This clearly advises against us- 
ing x 2 as a goodness-of-fit measure intended to 
rule out a particular model. 
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5.3. LSB rotation curve shapes 

The analysis procedures used by various au- 
thors induce significant differences in the rota- 
tion curves derived from the data, as illustrated 
by comparing the rotation curves of UGC 11557, 
and F568-3, the two galaxies common to the B01 
and S03 samples, and UGC 4325, the only galaxy 
common to the B02 and S03 datasets. (Figure 8). 
The rotation curve of UGC 11557 (Figure 8, upper 
left panel) presented by B01 extends out to only 
r ~ 6 kpc and rises with a nearly constant slope 
out to this radius. The best fit multi-parameter 
function has a scale radius r t = 14.1 kpc, much 
greater than the outermost radius of the obser- 
vations, r out cr = 6.2 kpc, and a low value of 
7 = 0.98. The S03 rotation curve (Figure 8, up- 
per right panel) for the same galaxy extends out 
to r ~ 10 kpc and appears to level off in the last 
two data points. The scale radius of the fit to this 
curve is r t = 6 kpc: well within the outermost 
radius of the observations r outC r = 10.4 kpc, and 
different by a factor of two from that of the B01 
curve. Although the value of the shape parameter 
7 is actually lower for the S03 (7 = 0.69), it is clear 
that additional data points in the flat part of the 
curve, or smaller error bars on the last few points, 
would force the fit to higher values of 7. This 
would also result in fits with a much lower value 
of the asymptotic velocity Vo, which presently has 
a value much greater than the maximum veloc- 
ity measured in either of the B01 or S03 rotation 
curves. 

In the case of F568-3, BOl's rotation curve (Fig- 
ure 8, middle left panel) rises monotonically out to 
the last point and as a result, the multi-parameter 
function fits quite well (x 2 rcd = 0.031). The S03 
rotation curve extends out to the same outermost 
radius as the B01 curve and both curves are con- 
sistent with one another within the observational 
errors. Due to the different methods of estimat- 
ing velocities, however, the B01 rotation curve is 
much smoother than that of S03. The S03 rotation 
curve (Figure 8, middle right panel) rises to a max- 
imum at r ~ 6 kpc then falls to V c ~ 0.9 V max at 
the outermost radius r out0 r = 11-2 kpc. The fit to 
this curve is further complicated by a deviant data 
point at r = 7 kpc with a velocity, V c ~ 0.7 V^n ax , 
much lower than that of neighbouring points. It 
is obvious that no simple fitting function can pro- 
vide a good fit to this rotation curve. The best fit 



multi-parameter function has a poor goodness-of- 
fit statistic Xrod = 1-5 and is characterized by an 
inordinately sharp turnover (7 = 25.7). We also 
note that the error bars in the B01 and S03 data 
are comparable in size (~ ±8 kms -1 ) except near 
the radius of the discrepant point in the S01 curve, 
where error bars in the B01 data are twice as large 
(~ ±16 kms™ 1 ). Clearly, the shape of the rota- 
tion curve as parameterized by fitting functions 
like the one given by eq. 7 is strongly influenced 
by the sensitivity of the observations, as well as by 
the method used to determine the rotation curve 
from the raw data. 

Perhaps the most interesting case is the rota- 
tion curve of UGC 4325. The B02 and S03 versions 
of the rotation curve appear qualitatively different 
and are not consistent with one another within the 
error bars as presented. The B02 rotation curve 
rises monotonically out the last measured point, 
whereas the S03 curve flattens off sharply to a 
well-defined asymptotic value. Fits to the B02 
and S03 rotation curves yield significantly different 
values for all three fitting parameters. The sharp 
turnover of the S03 curve results in a relatively 
large value of 7 = 3.67 compared to 7 = 1.38 for 
the fit to the B02 rotation curve. In addition, the 
asymptotic velocity Vo of the B02 fit is more than 
twice that of the S03 fit. Of even greater concern 
is the maximum velocity one might infer for UGC 
4325 from each version of the rotation curve. In 
the case of the B02 curve, the maximum velocity is 
undetermined but certainly appears to be greater 
than U(r outC r) — 122 kms -1 since the curve is still 
rising at the outermost data point. According to 
the S03 curve, however, the maximum velocity is 
robustly determined, equal to the asymptotic ve- 
locity of the fit, Vo — 94 kms -1 . The maximum 
velocity one infers for UGC 4352 therefore differs 
by at least 30% depending on which rotation curve 
is used. This may have important consequences 
for the calculation of scaling parameters such as 
the central densities of galaxies, a point we return 
to in §5.6. 

The spread in the best-fit parameters listed in 
the panels of Figure 8 is again a sobering reminder 
that a fair amount of non-trivial processing is in- 
volved when deriving rotation curves from raw 
data. Clearly, the disagreement between authors 
is somewhat worrying, and it limits the general 
applicability of conclusions inferred from individ- 
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ual galaxy studies. In spite of this, there appears 
to be broad statistical consistency between the ro- 
tation curve parameters derived by B01, B02 and 
S03. This is illustrated in the top panel of Fig- 
ure 9, where we show the distribution of best-fit 7 
values obtained for each sample. Each histogram 
in this figure is normalized to the total number of 
systems in each sample for ease of comparison. All 
three rotation curve datasets are broadly consis- 
tent with each other; in each case most (70%±5%) 
LSB rotation curves are characterized by a value 
of 7 < 2. These are typically gently-rising curves 
which turn over gradually as they approach the 
maximum asymptotic rotation speed, as shown, 
for example, by F574-1, and UGC 11454 in Fig- 
ure 6. A significant (~ 30%) number of 'outlier' 
galaxies with 7 ^ 2.5, however, populate the tail 
of the combined B01+B02+S03 distribution (see 
bottom panel of Figure 9). These are galaxies 
whose rotation curves feature a much "sharper" 
transition from the rising to flat part, as shown, 
for example, by F563-V2 and UGC 11748 in Fig- 
ure 6. 

5.4. Halo circular velocity profile shapes 

The bottom panel of Figure 9 shows the 7 dis- 
tribution obtained by fitting eq. 7 to the V c profile 
of all dwarf and galaxy-sized halos. In order to 
consider the various dynamical instances of a halo, 
we have included in the analysis about 20 differ- 
ent outputs for each system, spanning the redshift 
range 1 <j z < 0, giving a total of 266 halo profiles. 
We calculated the V c profile at each redshift in bins 
1 kpc (physical) in width for the galaxy halos and 
0.2 kpc (physical) in width for the dwarf halos, 
starting at the innermost reliably resolved radius 
r conv . The V c profiles were fit out to r = 20 kpc 
(physical) for the galaxy halos and r — 8 kpc 
(physical) for the dwarf halos, although we note 
that the fitting parameters are fairly insensitive 
to the region fitted provided that the curvature of 
the profile is well resolved. Since no formal error 
bars exist for the halo profiles (Poisson errors are 
negligible for the numbers of particles in these ha- 
los), we assign a uniform error of ±1 kms -1 to all 
points. 

The bottom panel of Figure 9 shows the dis- 
tribution of 7 values obtained for the halos, after 
convolution with the typical uncertainty in 7 de- 
rived from fits to the observed rotation curves (an 



asymmetric Gaussian with <r 7+ = 1.0, <t 7 _ = 0.5). 
The 7 distribution of all galaxy and dwarf halos, 
that of the dwarf halos only, and that of the com- 
bined B01, B02 and S03 sample are shown as the 
green, red and open histograms, respectively. The 
green and red histograms are normalized to the 
total number of halos, and the open histogram is 
normalized to the total number of galaxies in the 
combined observational dataset. As shown in Fig- 
ure 9, the convolved halo 7 distribution peaks at 
7 ~ 0.6, and has a dispersion of order ~ 0.6. For 
illustration, the three Gl V c profiles shown in the 
bottom panel of Figure 5 have 7 = 0.73, 0.65, and 
0.48 at z = 1.1, 0.48, and 0, respectively. There 
is no significant difference between the galaxy and 
dwarf halo 7 distributions. 

5.5. Comparison 

How does the distribution of rotation curve 
shapes compare with that of halo circular velocity 
profiles? The bottom panel of Figure 9 shows that, 
although there is significant overlap between the 
two distributions at values of 7 <j 2, most LSBs 
cluster around 7 ~ 1.2 compared to 7 ~ 0.6 for the 
halos. We note, however, that the contribution of 
a baryonic component has not been taken into ac- 
count in our analysis of the simulated halo V c pro- 
files. In order to investigate the effect of a baryonic 
disk on the shape of the V c profile, we construct 
an analytic mass model comprised of an NFW halo 
and an exponential disk. We use the prescription 
of Mo et al. (1998) to determine the scale length of 
the disk, Rd, as a function of the concentration, c, 
and spin parameter, A, of the NFW halo, and the 
mass, rrid, and angular momentum, jd, of the disk 
(expressed as fractions of the halo mass and angu- 
lar momentum, respectively). Fitting eq. 7 to the 
inner 20 kpc of the resulting V c profiles, we find 
that the best fit 7 value changes from 7 = 0.72 for 
an NFW halo with c = 10 and r 200 = 200 /i _1 kpc, 
to 7 = 0.92 (1.26) with the addition of a disk with 
m d = 3d = 0.05 (0.1) assuming A = 0.1, as ex- 
pected for LSB disks. 

We therefore conclude that the presence of the 
disk might be responsible for a significant sift in 
the peak of the 7 distribution. 

This suggests that the shape of most LSB ro- 
tation curves might actually be consistent with 
the circular velocity profiles of ACDM halos. The 
major difference in the two distributions comes 
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from considering the outliers; in particular, sys- 
tems with 7 ;> 2.5. Almost 1 in 3 LSBs is an 
outlier according to this definition, but such high 
values of 7 are rare amongst halos: fewer than 1 in 
40 halos have 7 > 2.5, and fewer than 1 in 100 have 
7 > 3. The "cusp vs. core" discrepancy alluded 
to above appears confined to less than a third of 
LSBs; those with sharp (7 ^ 2.5) turnovers in 
their rotation curves. 

We highlight the disagreement in Figure 10, 
where we show the V c profiles of all halos at z < 1 
alongside two rotation curves; one with 7 ~ 1 and 
another with a relatively high 7 ~ 5. In order 
to concentrate on the shape of the rotation curve 
rather than on the physical scaling parameters, 
profiles have been scaled in this plot to the ra- 
dius where the logarithmic slope of the fit to the 
rotation curve is equal to dlogV/dlogr = 0.3. 
We shall refer to this radius and its correspond- 
ing velocity as r .3 and V0.3, respectively. These 
parameters are easily retrieved from the fits with 
eq. 7, and are given by: r .3 = (7/3) _1 / 7 r t and 
Fo.3 = (10/7)- 1 /7^ . " 

The scaled halo profiles are described reason- 
ably well, on average, by the NFW profile, shown 
as the dashed curve in Figure 10, as is that of 
F571-8 (7 = 0.83). The shape of the rotation 
curve of F568-3 (7 = 5.4), on the other hand, is 
clearly inconsistent with the halo profiles. Better 
(albeit not perfect) fits to galaxies with 7 ;> 2.5 
may be obtained using the circular velocity curve 
of a system with a constant density core, as shown 
by the pseudo-isothermal 11 model indicated by the 
dot-dashed curve in Figure 10. 

5.6. The concentration of LSB halos 

The discussion of the preceding section focused 
on the shape of the rotation curves and halo V c 
profiles. We now turn our attention to the phys- 
ical parameters of the fits, in order to address 
claims that LSB galaxies are surrounded by ha- 
los of much lower concentration than expected in 
the ACDM scenario (McGaugh and de Blok 1998; 
dc Blok ct al. 2001b). We emphasize again that 
it is important to characterize both the observa- 



tional data and the simulations in a way that is 
as independent as possible from fitting formulae 
or extrapolation. Alam et al. (2002) recently pro- 
posed a simple and useful dimcnsionlcss measure 
of mass concentration that satisfies these criteria, 



V/2 



p(ry/2) 

Pcrit 



(8) 



11 The density profile of this widely used approximation to 
the non-singular isothermal sphere is given by pi so (r) = 
po/(l + ( r / r c) 2 , where r c is the core radius and po is the 
characteristic density of the core. 



Ay/2 measures the mean density contrast (rel- 
ative to the critical density for closure) within the 
radius at which the rotation speed drops to one 
half of its maximum value, V max . In practice, we 
estimate p(r) by 3 V c 2 (r) /AwGr 2 , a quantity that 
is easily measured both in galaxies with rotation 
curve data (and well defined Vmax) an d m simu- 
lated halos. 

The top panel of Figure 11 shows Ay/ 2 as a 
function of V ma x for all galaxies in the B01, B02, 
and S03 samples (open symbols), together with 
the galaxy halos in our sample (filled circles). 
Filled triangles and squares correspond to simu- 
lated dwarf galaxy and cluster halos, respectively 
(see § 2). 

The solid curve corresponds to the predictions 
of the Eke et al. (2001) halo concentration model 
for NFW halos in the ACDM cosmology we have 
assumed for our simulations. The dashed curves 
show the predictions of the Bullock et al. (2001) 
concentration model, together with the 1-cr halo- 
to-halo scatter predicted by their model. Figure 11 
shows that the simulations are in rough agreement 
with both models; the Bullock et al. (2001) model 
reproduces the simulations slightly better on the 
scales of dwarf halos, whilst the ENS model does 
better on the scale of clusters. 

Figure 11 also shows that, on average, the con- 
centration of ACDM is roughly comparable to 
those of LSBs in the samples we considered, al- 
though the scatter appears much larger than either 
obtained in the simulations or expected from the 
analytic model of Bullock et al. (2001). This result 
is reminiscent of our findings regarding the dis- 
tribution of rotation curve parameter 7 (see Fig- 
ure 9) : CDM halos appear consistent with the bulk 
of LSBs but are at odds with the most deviant sys- 
tems in the sample. 

In order to investigate the effect of observa- 
tional uncertainties on this parameter we exam- 
ine the values of Ay/2 

calculated for UGC 4325. 
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As discussed in §5.2, the maximum velocities mea- 
sured by B02 and S03 differ by 30% for this galaxy. 
The corresponding values of A v / 2 are shown as the 
open triangle (B02) and open square (S03) points 
in the top panel of Figure 11. Surprisingly, we find 
that the inferred halo central densities are very 
similar for the two versions of the rotation curve. 
Because Ay/2 depends only on the ratio between 
V max and the radius Ry/2, two curves with very 
different maximum velocities can give the same 
value of Ay/2 provided that the rising parts of the 
curves have the same slope. The halo central den- 
sities one infers from the B02 and S03 versions of 
the rotation curve for UGC 4325 are within 20% 
of one another despite a 30% difference in V" max 
between the two curves. This level of uncertainty 
in A v /2 is negligible considering that the halo-to- 
halo scatter in the simulations and model predic- 
tions spans a factor of 6 at a given V max . 

Although the effect of observational uncertain- 
ties may be relatively unimportant in the case 
of UGC 4325, we nonetheless attempt to reduce 
the scatter in the Ay/ 2 by culling from the sam- 
ple all galaxies whose rotation curves are still ris- 
ing at the outermost radius probed by the data. 
The bottom panel of Figure 11 shows the cen- 
tral densities inferred from rotation curves with 
rflog Vydlogr(r outor ) < 0.1. We find that this re- 
duced sample retains much of the original scatter, 
and many points lie both above and below the 
model predictions at a given value of V maK . 

Our results thus agree with those of Zentner 
and Bullock (2002), who have argued that the 
large number of low-concentration galaxies in LSB 
samples calls for substantial revision of the "con- 
cordance" ACDM scenario, such as tilted power 
spectra, running spectral index, or perhaps a lower 
erg. However, it appears unlikely that any such 
modification to the ACDM scenario would result 
in the larger variety of rotation curve shapes and 
concentrations apparently demanded by the LSB 
datasets. It is unclear at this point how to reduce 
the disagreement, but any resolution to the puzzle 
must explain why so many LSBs are actually in 
good agreement with ACDM halos and why the 
disagreement is confined to a minority of systems. 
The possibility remains that some complex astro- 
physical process not yet considered in the models 
might actually be behind the discrepancy and that 
no radical modification to the ACDM paradigm is 



called for. 

6. Conclusions 

We present results from a set of high resolu- 
tion cosmological simulations of dark matter halos 
formed in a ACDM cosmogony. Seven Milky Way- 
sized galaxy halos were simulated at various mass, 
time, and spatial resolutions, enabling us to inves- 
tigate the convergence properties of cosmological 
N-body simulations. We have examined the inter- 
nal structure of the highest resolution realization 
of each halo, with particular emphasis on the log- 
arithmic slope of the inner density profile. Finally, 
we have compared the circular velocity profiles of 
these halos with the rotation curves of a large sam- 
ple of dwarf and LSB galaxies. 

Our main conclusions may be summarized as 
follows. 

• The convergence criteria proposed by Power 
et al. (2003) are robust, and provide a con- 
servative estimate of the minimum radius 
at which the mass profile of simulated ha- 
los can be reliably predicted. According to 
these criteria, the highest resolution galaxy 
halos we have simulated (which contain 2-4 
million particles within the virial radius) are 
reliably resolved down to r conv ~ 1 /i _1 kpc. 

• The slope of ACDM halo density profiles 
becomes progressively shallow all the way 
down to the minimum reliably resolved ra- 
dius, without sign of converging to a a well- 
defined power law near the centre. 

• In general, the slope changes with radius 
more gradually than predicted by the NFW 
formula, which leads some halos to be bet- 
ter described by profiles with steeper cusps, 
such as the modification to the NFW for- 
mula proposed by Moore et al (1999). There 
is, however, significant variation from halo to 
halo in the radial dependence of the slope. 
Some systems are better fit by the NFW 
profile, and others by the Moore et al for- 
mula. At r conv , however, the density pro- 
files are significantly shallower than r -15 , 
the asymptotic value advocated by Moore 
et al. (1999). 
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• A comparison of the circular velocity pro- 
files of ACDM halos with rotation curves of 
dwarf and LSB galaxies indicates that the 
shapes of most (~ 70%) LSB rotation curves 
are consistent with those of simulated halo 
V c profiles. The remaining ~ 30% feature 
a sharp turnover from the rising to the flat 
part of the curve which is not consistent with 
the structure of simulated halos. 

• The concentration of dwarf and galaxy- 
sized halos simulated in the "concordance" 
ACDM cosmogony is in reasonable agree- 
ment with the characteristic central density 
inferred for most LSB galaxies from rotation 
curve data. However, as with the rotation 
curve shape, a significant number of LSBs 
have concentrations well below (and above) 
the expected range. 

We conclude that the inner structure of ACDM 
halos is not manifestly inconsistent with the rota- 
tion curves of LSB galaxies, although they seem 
unable to reproduce the full variety of LSB ro- 
tation curve shapes and normalizations. This dis- 
crepancy may signal the need to revise some of the 
basic tenets of the ACDM scenario, but it might 
also taken to imply that the relation between gas 
rotation speeds and spherically-averaged halo cir- 
cular velocities is more complex than assumed in 
simple analyses such as the one presented here. 

CDM halos, for example, are known to be tri- 
axial, which may lead gaseous disks to deviate sys- 
tematically and significantly from simple coplanar 
circular orbits. Work is in progress to try and 
determine whether such asymmetries in the po- 
tential are able to account quantitatively for the 
"cusp vs core" discrepancy. Until this or other 
plausible astrophysical origin for the discrepancy 
is identified, the observed variety of LSB rotation 
curves and concentrations will remain a challenge 
that the theory must meet. 
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Table 1 

Numerical and Physical Properties of Simulated Halos 
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Table 1 — Continued 



Label 
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Code 
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Table 2 

Properties of Rotation Curves and Fit Parameters 
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Table 2 — Continued 
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Fig. 1. — The structure of halo G3 at different mass, time, and spatial resolutions. The value of the 
softening parameter e is indicated by arrows in the three panels. The number of timesteps and particles 
are listed in Table 1. Runs with 32 3 , 64 3 , 128 3 and 256 3 high-resolution particles are shown in blue, 
green, red, and black, respectively. Dotted curves show an NFW profile with concentration c = 5.3 and 
^200 = 143.4 /i~ 1 kpc. Upper panel: Circular orbital timescale i C i rc versus radius; radii at which the circular 
orbital timescale is less than 15 N At ' , indicated by the dashed lines for each simulation, are unresolved 
due to insufficient time resolution. Middle panel: Mean radial acceleration profile V c (r) 2 /r; untrustworthy 
radii are those corresponding to accelerations greater than the limiting acceleration imposed by the softening 
a e ~ 0.5 V 2 oo/e, shown by the dashed lines. Lower panel: Collisional relaxation time i re iax versus radius; 
convergence requires £ rc iax icirc(T20o)- Vertical dotted lines indicate the radius, r conv , beyond which this 
condition (the strictest of the three) is satisfied. 
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Fig. 2. — Upper left panel: Density profiles of halo G3 at four different levels of mass resolution, plotted 
down to radii containing 50 particles. Dashed and dotted curves show best fit Moore et al. and NFW profiles, 
respectively (see text for fitting details). Vertical dotted lines indicate the minimum converged radius, r GO nv, 
for each run. Upper right panel: Same density profiles plotted only for converged radii. The discrepancy 
between lower resolution runs and the highest resolution simulation at small radii is no longer apparent when 
only reliably resolved radii are considered. Lower panels are as the upper ones, but for the circular velocity 
profiles. 
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Fig. 3. — Radius where the circular velocity profile of lower resolutions runs starts to deviate from that of the 
highest resolution run by more than 10%, r w y oVC , plotted against the minimum converged radius predicted 
by the P03 convergence criteria, r conv . In all cases r conv < r w % vci validating the P03 criteria as conservative 
estimators of the converged region. 
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Fig. 4. — Upper left panel: Logarithmic slope of the smoothed density profile of halos simulated at our highest 
resolution N s \, ox = 256 3 , plotted for r > r conv . Curves are scaled horizontally to the radius r_ 2 , where the 
slope takes the isothermal value, dlog/d/(ilogr(r_2) = —2. NFW and Moore et al. profiles are shown as solid 
and dashed curves, respectively. The logarithmic slope increases monotonically with decreasing radius and 
there is no obvious convergence to a particular asymptotic value of the central slope. Upper right panel: Same 
as upper-left but for the SCDM Virgo cluster of Ghigna et al. (2000) and SCDM Milky Way (MW) and M31 
galaxy halos of Moore et al. (1999). The profiles of the two galaxy-sized halos appear to be consistent with 
those of our halos. The logarithmic slope of the cluster halo appears to be slightly steeper than the others, 
fluctuating about a value of —1.4 at the innermost resolved point. Lower left panel: Halo density profiles 
scaled horizontally to radius r_2 and vertically to the corresponding density at that radius p_2 = p{r-2)- 
Lower right Panel: Same as lower- left but for the Ghigna et al. (2000) and Moore et al. (1999) halos. 
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Fig. 5. — Top panel: Mass within r = 2, 10, and 20 kpc (physical) and r2oo for halo Gl/256 3 as a function 
of redshift (age of Universe in Gyr) on bottom (top) axis. The mass within 20 kpc undergoes significant 
fluctuations in response to a merger at z ~ 0.7. Bottom panel: The inner circular velocity profile before 
(z — 1.10), during (z = 0.48), and after (z = 0) the merger. The shape of the V c profile is noticeably altered 
by the effects of the infalling substructure. 94 




5 10 15 200 5 10 15 200 5 10 15 20 



r (kpc) r (kpc) r (kpc) 



Fig. 6. — LSB rotation curves from the datasets of de Blok et al. (2001b) (B01), de Blok and Bosma (2002) 
(B02) and Swaters et al. (2003) (S03). Solid curves show best fits using the Courteau (1997) fitting formula 
given by eq. 7. Best fits with 7 *~ 1 correspond to rotation curves that rise and turn over gently as a function 
of radius. Curves with 7 > 2 feature a much sharper transition from the rising to the flat part of the curve. 
See text for full discussion. 
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Fig. 7. — Distribution of reduced x 2 values corresponding to best fits using the multi-parameter fitting 
formula given by eq. 7 for the three rotation curve datasets (B01, B02, S03). Note that the B01 and B02 
Xied distributions peak at unrealistically low values, signalling the presence of significant correlation between 
neighbouring points in the rotation curves of these samples. The S03 % 2 distribution peaks at a higher value, 
emphasizing the various assumptions of different authors when deriving rotation curves from raw data (see 
also Figure 8). 
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Fig. 8.— Rotation curves of UGC 11557, F568-3, and UGC 4325 derived by either de Blok et al. (2001b); de 
Blok and Bosnia (2002) (BOI and B02, left panels) or Swaters et al. (2003) (S03, right panels). Solid curves 
show fits to the rotation curves using eq. 7. Values of the three fitting parameters 7, rt, and Vq are listed in 
each panel along with the corresponding reduced chi-squared, xj? cd . This figure illustrates the effect on the 
derived rotation curves due to different assumptions and methodology adopted by various authors. In the 
case of UGC 4325, for example, the rotation curve derived by B02 continues to rise out to the last measured 
point at V ~ 120 kms -1 , whereas the the S03 curve flattens off at V ~ 90 kms -1 . 
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Fig. 9. — Top panel: Distribution of best-fit 7 values obtained for galaxies in the samples of de Blok 
et al. (2001b) (B01), de Blok and Bosnia (2002) (B02), and Swaters et al. (2003) (S03). Note that most 
galaxies in all three samples cluster around 7 = 1.2, but that there are a significant number of outliers 
with 7 J> 2.5. These define a population of galaxies that seems to be distinct from the bulk of the sample. 
Bottom panel: Combined observational sample compared with the halo 7 distribution after convolution with 
an error distribution similar to that corresponding to observed rotation curve fits. Arrows show the change 
in 7 caused by the addition of an exponential disk to an NFW halo with concentration c = 10. The mass 
of the disk, rrid, is given in units of the halo mass. Its exponential scale length is computed assuming that 
the spin parameter of the halo is A = 0.1, and that the specific angular momentum of the disk is the same 
as that of the halo (Mo et al. 1998). The magnitude of the correction (shown by horizontal arrows) suggests 
that the 7 distribution of halos might actually be consistent with that of the bulk of galaxies, but apparently 
fails to account for the 7 J> 2.5 tail in the rotation curve distribution. See text for further details. 
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Fig. 10. — Rotation curves of two galaxies, one with 7 = 0.83 (f571-8, filled circles), representative of the 
bulk of the population, and another one with 7 = 5.4 (f568-3, filled triangles) representative of the tail of 
the distribution. The curves have been scaled to the radius, ro.3, (and corresponding velocity, V0.3) where 
the logarithmic slope of the fit to the rotation curve is d\ogV/d\ogr = 0.3. Galaxy (dwarf) halos are shown 
by solid (dotted) lines, and they appear to follow closely the NFW profile (dashed line). This figure shows 
clearly that halo V c profiles are inconsistent with rotation curves with high 7, such as f568-3. The dot-dashed 
curve shows that f568-3 is somewhat better matched by a non-singular isothermal sphere model. 
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Fig. 11. — Alam et al. (2002) halo concentration parameter, Ay/2) inferred from rotation curves of LSB 
galaxies in the B01, B02, and S03 samples (open symbols), compared to those measured for simulated dwarf- 
(logVmax — 1.6), galaxy- (log Vmax — 2.3), and cluster-sized (log Vmax — 3.2) halos (solid symbols). Points 
corresponding to Ay/ 2 values for UGC 4325 calculated from B01 and S03 rotation curves are marked with 
an x. Solid line shows the prediction of the Eke et al. (2001) concentration model for NFW halos in a 
ACDM cosmogony (f2 = 0.3, = 0.7, and h — 0.65, a§ = 0.9). Dashed lines show the prediction (and 1 a 
scatter) corresponding to the Bullock et al. (2001) concentration model. Top panel: All galaxies in combined 
B01, B02, and S03 dataset. Bottom panel: Only galaxies whose rotation curves are not rising steeply at the 
outermost point (dlog V/dlogr(7" utcr) < 0.1). 
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